!crop***********************************************************************
      subroutine crop(jd,hr,tair,pgros,ww,stage,hsum,tsenlai, &
          tlai,pg,rm,dw,newlai,emergedate,droughtflag,sstart,cutcnt, &
          lfwt,stwt,agb,grwt,harvest,rtwt)
!
!     Samuel Adiku's Crop Model
!--------------------------------------------------------------------------
      IMPLICIT NONE


      integer ::  nstag

      parameter (nstag=5)

      common/crops/sowdate,cutday1,cutday2,cutday3,grassflag,senstage
      common/cropss/hs,conveff,minbiom,cutlai,startlai,partl,partst, &
          partrt,partgr,SLA,senfac,efv,laytck,maxrtdep,tb

      double precision tb,hs(nstag),conveff,partl(nstag),partst(nstag), &
          partrt(nstag),partgr(nstag),SLA,senfac,efv,laytck, &
          maxrtdep,minbiom,cutlai,startlai
      integer sowdate,emergedate,cutday1,cutday2,cutday3, &
          grassflag,senstage,droughtflag,sstart,send,cutcnt

      integer jd,dap,nrlays,stage
      integer harvestday,harvest,harvestcnt
      double precision hr,tair,pgros,pg,rm,dw
      double precision dtt,dlfwt,dstwt,drtwt, &
               dgrwt,ndlai,hsum, &
               crtdep,tlai,tsenlai,dsenlai
      real ww,newlai,lfwt,stwt,agb,grwt,rtwt

      save dap,harvestday,harvestcnt    !rtwt,grwt????

!2dim hsum,dap,rtwt,grwt,tlai

      print*, 'subroutine crop'


      if (jd.eq.1.and.hr.eq.0.) then
         harvestday=1000
         harvestcnt=0
      endif

      call heats(jd,hr,tair,tb,sowdate,dtt)
      call crstage(hsum,dap,hs,stage)

      if (stage.gt.2.and.harvestcnt.gt.0) stage=2

      call growth(tair,conveff,pgros,pg,rm,dw,ww)
!wet
      send=395-sstart
      if (harvestday.lt.400) then
         if (droughtflag.eq.0.and.jd.gt.send) dw=0.
      endif

      call plpart(partl(stage),partst(stage), &
                 partrt(stage),partgr(stage), &
                 dw,dlfwt,dstwt,drtwt,dgrwt)

!     calculates the increment in leaf area
      ndlai=dlfwt*SLA/10000.    ! 10000 converts from cm2 to m2

      dsenlai=0.d0
!wet
      if (stage.ge.senstage.and.harvestcnt.eq.0) &
          dsenlai=senfac*tlai

      if (droughtflag.eq.0.and.jd.gt.send.and. &
          harvestcnt.gt.0) then
         dsenlai=newlai/50./24.
         dlfwt=dlfwt-dsenlai/SLA*10000. ! dlfwt and dstwt should be 0 here, because dw=0
         dstwt=dstwt-dsenlai/SLA*10000.
      endif

      call rtdep(dap,efv,crtdep,nrlays,laytck,maxrtdep)
      call rtgrowth(stage,drtwt,nrlays)

!rtsen-------------------------------------------------------------------------
      if (stage.ge.senstage.and.harvestcnt.eq.0) &
          drtwt=drtwt-dsenlai/SLA*10000.
!rtsen-------------------------------------------------------------------------

      hsum=hsum+dtt

      ww=ww+dw          ! total or cumulative biomass

      lfwt=lfwt+dlfwt   !gC m-2
!wet
      if (lfwt.lt.0.) lfwt=0.
      stwt=stwt+dstwt

      rtwt=rtwt+drtwt

      grwt=grwt+dgrwt

      if (ndlai.le.0.) ndlai=0.

      tlai=tlai+ndlai
      tsenlai=tsenlai+dsenlai

      if (tsenlai.gt.tlai) then
         tsenlai=tlai
         if (harvestday.gt.500) then
            harvestday=jd+14    !2 weeks after LAI=0
            harvest=harvestday
         endif
      endif

!     to force harvest
      if (jd.eq.300.and.harvestday.gt.500) then
         harvestday=jd+2
         harvest=harvestday
      endif

      newlai=tlai-tsenlai

      if (newlai.gt.1..and.harvestcnt.gt.0) then
         newlai=1.
         tlai=newlai
      endif

      if (jd.eq.emergedate) tlai=startlai
      if (jd.ge.sowdate) dap=jd-sowdate

      if (jd.eq.harvestday) then !harvest
         stage=0
         ww=0.
         lfwt=0.
         stwt=0.
         grwt=0.
         rtwt=0.
         tlai=0.1
         newlai=tlai
         hsum=hs(1)
         tsenlai=0.
         harvestcnt=1
      endif

      return
      end










!crop***********************************************************************
      subroutine crstage(hsum,dap,hs,stage)
!
!     determines stage of the crop
!---------------------------------------------------------------------------
      implicit none

      integer, parameter :: nstag =5

      integer stage,dap
      double precision hsum,hs(nstag)

      if (dap.ge.1) then
         if (hsum.le.hs(1)) stage=1
         if (hsum.gt.hs(1).and.hsum.le.hs(2)) stage=2
         if (hsum.gt.hs(2).and.hsum.le.hs(3)) stage=3
         if (hsum.gt.hs(3).and.hsum.le.hs(4)) stage=4
         if (hsum.gt.hs(4).and.hsum.le.hs(5)) stage=5
      endif

      return
      end






!crop***********************************************************************
      subroutine growth(tair,conveff,pgros,pg,rm,dw,ww)
!
!     calculates dry weight (dw) from pgros
!--------------------------------------------------------------------------
      implicit none

      double precision pgros,dw,Rm,km,tair,pg,conveff
      real ww

      km=0.0006             ! (see Jones 1991)
!    Ec=conveff=0.70       ! conversion efficiency from photosynthesis to dry matter
!     convert mmol m-2 h-1 -> gC m-2 h-1
!     gpp=pgros/1000.*12.

!     convert umol CO2 m-2s-1 to g CH2O m-2h-1
      pg=pgros*0.108        ! 0.108 is conversion factor

!     Now calculate growth and maintenance respiration (Ref. Jones, 1991)
!     First maintenace respiration rm, as a function of air temperature
      rm=km*exp(0.0693*(tair-25.))

!     Now determine crop dry weight increment
      dw=conveff*(pg-rm*ww)   ! Ec accounts for growth respiration as well

      return
      end






!crop***********************************************************************
      subroutine plpart(pl,ps,pr,ppg,dw,dlfwt,dstwt,drtwt,dgrwt)

!     determines partitioning coefficients
!--------------------------------------------------------------------------
      implicit none

      double precision pl,ps,pr,ppg,dw,dlfwt,dstwt,drtwt,dgrwt

      if (dw.ge.0.) then
         dlfwt=pl*dw
         dstwt=ps*dw
         drtwt=pr*dw
         dgrwt=ppg*dw
      else
         dlfwt=0.
         dstwt=0.
         drtwt=0.
         dgrwt=0.
      endif

      return
      end






!crop***********************************************************************
      subroutine rtdep(dap,efv,crtdep,nrlays,laytck,maxrtdep)

!     This subroutine calculates the current root depth
!     Need to improve the concept here or find actual efv values
!--------------------------------------------------------------------------
      implicit none

      integer dap,nrlays
      double precision efv,maxrtdep,crtdep,laytck,frrtd,diff

!     efv=1.5                ! cm/d
!     efv=2.0
!     elseif (mk.eq.2) then  ! maize
!     efv=2.
!     maxrtdep(mk)=200.
!     else
!     maxrtdep(mk)=100.      ! for grassland
!     efv(mk)=1.
!     endif

      crtdep=efv*dap

      if (crtdep.ge.maxrtdep) crtdep=maxrtdep

!     gy: better use modulo operator
      frrtd=crtdep/laytck
      nrlays=int(frrtd)
      diff=frrtd-nrlays
      if (diff.lt.0.5) then
         diff=0.
      else
         diff=1.
      endif

      nrlays=nrlays+diff
      if (nrlays.le.1) nrlays=1

      return
      end







!crop***********************************************************************
      subroutine rtgrowth(stage,drtwt,nrlays)
!
!     This subroutine determines the root length density in each soil layer
!--------------------------------------------------------------------------
      implicit none

      integer nrlays,i,stage
      double precision drtwt,initrld,rld(15),rtcon, maxrld,wf(15), &
          today_rld,delt_rld(15),tdrld,drld(15),envf


      initrld=0.5
      rtcon=6500.               ! cm/g
      maxrld=4.                 ! 4 cm cm-3
      envf=1.                   ! Environmental factor set to unity. This must be linked to water content

!     drwt is in g/m2 -> change to g/cm2 then multiply with rtcon and finnaly
!     convert to cm/cm3
      today_rld=rtcon*drtwt/100000.   ! converted to cm/cm3

      if (stage.ge.2) then
         tdrld=0.
         do i=1,nrlays
            if (rld(i).eq.0) rld(i)=initrld            ! .eq.0 fuck
            if (rld(i).ge.maxrld) rld(i)=maxrld
            delt_rld(i)=rld(i)*(1.-rld(i)/maxrld)*envf ! combined carrying capacity and environmental stress
            tdrld=tdrld+delt_rld(i)
         enddo
         do i=1,nrlays
            if (tdrld.gt.0.0) wf(i)=delt_rld(i)/tdrld
            drld(i)=wf(i)*today_rld
            rld(i)=rld(i)+drld(i)
         enddo
      endif

      return
      end




